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ABSTRACT 

In this paper, the third in a series illustrating the power of generalized linear models 
(GLMs) for the astronomical community, we elucidate the potential of the class of 
GLMs which handles count data. The size of a galaxy’s globular cluster population 
(Ngc) is a prolonged puzzle in the astronomical literature. It falls in the category 
of count data analysis, yet it is usually modelled as if it were a continuous response 
variable. We have developed a Bayesian negative binomial regression model to study 
the connection between Nqq and the following galaxy properties: central black hole 
mass, dynamical bulge mass, bulge velocity dispersion, and absolute visual magnitude. 
The methodology introduced herein naturally accounts for heteroscedasticity, intrinsic 
scatter, errors in measurements in both axes (either discrete or continuous), and allows 
modelling the population of globular clusters on their natural scale as a non-negative 
integer variable. Prediction intervals of 99 per cent around the trend for expected Nqc 
comfortably envelope the data, notably including the Milky Way, which has hitherto 
been considered a problematic outlier. Finally, we demonstrate how random intercept 
models can incorporate information of each particular galaxy morphological type. 
Bayesian variable selection methodology allows for automatically identifying galaxy 
types with different productions of GCs, suggesting that on average SO galaxies have 
a GC population 35 per cent smaller than other types with similar brightness. 

Key words: methods: statistical, data analysis-galaxies: globular clusters 


1 INTRODUCTION 


The current era of astronomy marks the transition from a 
data-deprived field to a data-driven science, for which sta¬ 
tistical methods play a central role. An efficacious data ex¬ 
ploration requires astronomers to go beyond the traditional 
Gaussian-based models which are ubiquitous in the field. 
Gaussian distributional assumptions fail to hold when the 
data to be modelled come from exponential family distri¬ 


butions other than the Normal/GaussiarQ (|Hardin &; Hilbe 


2012 Hilbe 20141. For non-Gaussian regression problems 


there exist powerful solutions already widely used in medi¬ 
cal research (e.g., |Lind scy 199 9]) , finance (e.g., |Jong fe Heller| 
20081, healthcare (e.g., Griswold et al.|[2004 1 and biostatis¬ 


tics (e.g., Marschner & G illett| |2012[ ), but vastly under¬ 
utilized to-date in astronomy. These solutions are known as 
generalized linear models (GLMs). Despite the ubiquitous 


1 The exponential family comprises a set of distributions ranging 

from both continuous and discrete random variables (e.g., Gaus¬ 

sian, Poisson, Bernoulli, Gamma, etc.) 
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implementation of GLMs in general statistical applications, 
there have been only a handful of astronomical studies ap- 


plying GLM techniques such as logistic regression (e.g. Rai- 

choor & Andreon 2012| 2014| Lansbury et al. 2014 ; De Souza 

et al. |2015), Poisson regression (e.g 

. Andreon & Hurn|2010), 

gamma regression (Elliott et al. 

10151 and negative bino- 

mial (NB) regression (Ata et al. 

2015). The methodology 


discussed herein focuses on Bayesian count response mod¬ 
els (Poisson and NB), suited to handle discrete, count-based 
data sets applied to a catalogue of globular clusters (GCs). 

Globular clusters are among the oldest stellar systems 
in the Universe (formed at 2 > 2, Kruijssen|2014 1 , are perva¬ 
sive in nearby massive galaxies, |B rodie fc Strader] 2006) and 
can be found in massive galaxy clusters not necessarily asso¬ 
ciated to one of its galaxies (e.g., Durrell et al. |2014| . Hence, 
understanding their properties is of utmost importance for 
drawing a complete picture of galaxy evolution. The past 
few decades have seen considerable interest in the apparent 
correlation between the mass of the black hole at the centre 
of a galaxy, Mbh, and the velocity dispersion of the central 
stellar bulge, a (e.g., [Gcbhardt et al. |2000| . As part of the 
process of understanding the nature and origin of the so- 
called Mbh-u relation, astronomers have investigated links 
between other properties of the host galaxy. In particular, 
the correlation between the size of globular cluster popula¬ 
tions, Ngc, and Mbh is tight, possibly more so than the 
Mbh-u relation, and may reflect an underlying connection 
to the bulge mass, binding energy, host galaxy stellar mass 
and total luminosity (Burkert & Tremaine||2010| Harris & 
Harris||2011| |Snyder et al.||2011| |Rhode||2012| |Harris et al. 
20T3[ 2014JT This may go some way to explaining the huge 
range in scales of the regions involved. One notorious outlier 
is our own Milky Way galaxy, for which there are far too 
many globular clusters given the mass of its central super- 
massive black hole, despite the fact that both are accurately 
measured. Nevertheless, the otherwise small scatter found in 
such relations deserves a closer look since it cannot be easily 
explained by simple scaling relation arguments. 

The connection between Nqc and the global proper¬ 
ties of their host galaxies is an extant astronomical puzzle 
involving count models, but is treated as a continuous one. 
Such correlation studies are commonly based on taking pairs 
of parameters (x,y) in log-log space and searching for solu¬ 
tions in the normal form y = a + /3x, despite the fact that 
this regression technique assumes continuous variables and a 
Gaussian error distribution, e.g. ^-minimisation (Tremaine 
et al. |2002 1. 


Our method surpasses the previous x -minimisation ap¬ 
proach in several ways. The most obvious being the ability 
to handle count data without the need of logarithmic trans¬ 
formations of a discrete variable. Hence, we can take into ac¬ 
count the cases with zero counts, instead of removing them 
to accommodate the logarithm transformation, or adding 
an arbitrary data shift in the form log(x+e), with e com¬ 
monly taken as unity. Our method naturally handles errors 
in variables in both the x and y axes accommodating the 
heteroscedasticity of the errors in As a further anal- 


2 Heteroscedastic error structures may remain even after trans¬ 
formation, thus violating the Gaussian assumption of homogene¬ 
ity of error variance. 


ysis, we introduce one of the most important extensions of 
GLMs known as generalized linear mixed models (GLMMs). 
This is done to include in the model information about each 
galaxy morphological type, allowing discrimination among 
classes of objects requiring additional adjustments in their 
regression coefficients. 

The outline of this paper is as follows. In section [2] we 
provide a brief introduction of generalized linear models in 
the context of exponential family distributions. An overview 
of count data along with Poisson and NB GLMs are pre¬ 
sented in section [3] The dataset used in our analysis is sum¬ 
marized in section [4] In section [5] we discuss the necessary 
steps to build our Bayesian model. In section [6] we discuss 
GLMMs in the context of random intercepts models. Finally 
in section [7] we present our conclusions. 


2 GENERALIZED LINEAR MODELS 


Classical response-with-covariates models, that is, general 
(not generalized) linear models, assume that the response 
variable and the residual errors, following a normal distri¬ 
bution, are linear in the model parameters and have constant 
variance. This allows model parameter estimation with or¬ 
dinary least squares (OLS) methods. As described above, 
many data sets have response variables that violate one or 
more of these assumptions. While remedial measures such as 
transformations on the response variable or the covariates 
may be applied, these measures may fall short of satisfy¬ 
ing the OLS requirements. For data sets for which classical 
models are ill-suited, the extended class of models, GLMs, 
are used with model parameters often estimated using maxi¬ 
mum likelihood methods (for a brief overview of GLMs in an 
astronomical context, see e.g., De Souza et al.|[20T5| |Elliott| 
et al.|2015l. 


Nelder & Wedderburn (1972) introduced an unifica¬ 


tion of models characterised by being linear on the sys¬ 
tematic component (model predictors). For example logis¬ 
tic and probit analysis for binomial variates, contingency 
tables for multinomial variates, and regression for Poisson- 
and gamma- distributed variates, each a form of the GLM. 
The random response variable, Y t . i = 1, 2,..., n, may be 
represented as 


g{yi) = yt, 

rn = xf f3 = /3 0 + (hxi H-b /3pX p 


(1) 


In equation / denotes a response variable dis¬ 
tribution from the exponential family (EF), fM is the re¬ 
sponse variable mean, (j> is the EF dispersion parameter 
in the dispersion function a(-), V(im) is the response vari¬ 
able variance function, r/i is the linear predictor, the xf = 
{xn, Xa, ..., Xi p } T is the vector of explanatory variables 
(covariates or predictors), /3 = {/3i, 02 , ■ ■ ., flp} is the vec¬ 
tor of covariates coefficients, and g(-) is the link function, 
which connects the mean to the predictor. If V(fii) is a con¬ 
stant for all yi, then the mean and variance of the response 
are independent, which allows using a Gaussian response 
variable. If the response is Gaussian, then <?(/.i) = y. The 
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general form of the GLM thus allows Gaussian family, N, 
linear regression as a subset, taking the form: 


'-' 

Yi ~ A f(iM, a 2 ), 

!M = /3o + P 1 X 1 H-+ flpXp. 


( 2 ) 


The subset of GLMs for count data are the Poisson regres¬ 
sion models and the several incarnations of the NB regres¬ 
sions. Poisson regression models assume the count response 
variable follows a Poisson probability distribution function. 
Similarly, the NB regression models assume the count re¬ 
sponse variable follows a NB probability distribution func¬ 
tion. Descriptions of the Poisson and NB models follow. 


3 MODELLING COUNT DATA 


Astronomical quantities can be measured on different scales: 
nominal (e.g., classes of objects: Type Ia/II supernovae, el¬ 
liptical/spiral galaxies); ordinal, (e.g., ordering planets ac¬ 
cording to their size or distance to the star); and metric 
(e.g., galaxy mass, stellar temperature). Observations that 
have only right-skewed, non-negative integer values belong 
to a subclass of the metric scale known as count data. Dis¬ 
tances between counts are meaningful, hence the counts are 
metric, but they are not continuous and must be treated 
as such. Astronomical count data are often log-transformed 
to satisfy Gaussian parametric test assumptions rather than 
modelled on the basis of a count distribution. Despite the 
fact that GLMs are better suited to describe count data, a 
log-transformation of counts has the additional problem of 
dealing with zeros as observations. With just one observa¬ 
tion with value zero, the entire data set needs to be shifted 
by adding an arbitrary value before transformation. It is well 
known that such transformations perform poorly, leading to 
bias in the estimated parameters (O’H ara fc Kotzep OlO). 

We begin our discussion of regression models for count 
data with the subset of GLMs known as Poisson regression. 
A common condition accompanying count data is overdis- 
person, it occurs when the variance exceeds the mean. This 
condition in Poisson regression suggests that remedial mea¬ 
sures, such as the use of NB regression, may be appropriate. 


3.1 Poisson Regression 


Poisson regression was the first model specifically used to 
deal with count data and still stands as basis for many types 
of analyses. It assumes a discrete response described by a 
single parameter distribution which represents the mean or 
rate, /x; i.e., the expected number of times an event occurs 
within a fixed time-interval. Another important feature is 
the assumption of equidispersion which implies the equality 
of mean and variance, and c an b e quantified by the Pearson 
X 2 dispersion statistic (see 


3.2 


). The Poisson distribution 


function is typically displayed as 


f(y;n) = 


y'- 


where the mean and variance are given by 
Mean = /x, Variance = yu, 


(3) 


(4) 
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representing a particular case of equation 0 with V (/x) = /x 
and a(4 >) = 1. Thus, a regression equation derived from 
equation 0 may be used as a GLM for a count response, 
y. The usual link function, <?(/x), is the natural log func¬ 
tion such that /x = e 11 (see e.g., |Hardin fe Hi lbc 2012). It is 
worth noting that GLMs are not simple log transforms of 
the response variable, but rather, the expected counts from 
a Poisson regression is an exponentiated linear function of r/, 
thereby keeping the response variable on its original scale. 
Often, count data do not enjoy the Poisson assumption of 
equidispersion resulting in a Poisson dispersion statistic (see 
section [3(2)1 with a value greater than one. 


3.2 Overdispersion 


Overdispersion in Poisson models occurs when the response 
variance is greater than the mean. It may arise when there 
are violations in the distributional assumptions of the data 
such as when the data are clustered, thereby violating the 
likelihood requirement of the independence of observations. 
Overdispersion may cause standard errors of the estimates 
to be deflated or underestimated, i.e. a variable may ap¬ 
pear to be a significant predictor when it is in fact not. A 
key approach for checking overdispersion is by means of the 
dispersion statistic, V , 


77 = 


N-Np’ 


(5) 


where N is the number of observations and N p is the number 
of parameters in the model. Then N — N p represents the 
residual degrees of freedom. For a Poisson GLM, the Pearson 
X 2 value is 


X 


2 


E 


(Yj - 


( 6 ) 


where Y, represents the observed values, and Hi is the mean 
and variance of Y- Poisson overdispersion occurs when the 
variation in the data exceeds the expected variability based 
on the Poisson distribution, resulting in 77 being greater than 
1. Small amounts of overdispersion are of little concern; a 
rule of thumb is: if 77 > 1.25, then a correction may be 
warranted | Hilbe|[20T4 1. 

If overdispersion is observed, then there are several cor¬ 
rective measures in common practice. Options are adjusting 
the standard errors by scaling, applying sandwich or robust 
standard errors, or bootstrapping standard errors for the 
model. However, only the standard errors will be adjusted 
and not the regression coefficients, /3, which often can be 
affected by overdispersion as well (e.g., |ffiIbel|20TTl ). This 
paper examines the efficacy of using Bayesian estimation 
methods on a more general discrete distribution known as 
the NB. The NB distribution contains a second parameter 
called the dispersion or heterogeneity parameter which is 
used to accommodate Poisson overdispersion as described 
below. 


3.3 Negative Binomial Regression 

The NB distribution has long been recognized as a full mem¬ 
ber of the exponential family, originally representing the 
probability of observing y failures before the rth success in 
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Table 1. Main assumptions of each regression model family. 



Normal 

Log-normal 

Poisson 

Negative binomial 

Response variable 

Real 

Positive 

Non-negative integer 

Non-negative integer 

Null values 

✓ 

X 

/ 

/ 

Sample variance 

Homoscedastic 

Homoscedastic 

Heteroscedastic 

Heteroscedastic 

Overdispersion 

X 

X 

X 

/ 


a series of Bernoulli trials. It can also be formulated as a 
Poisson model with gamma heterogeneity (Hilbe 2011). The 
NB model, as a Poisson-gamma mixture model, is appropri¬ 
ate to use when the overdispersion in an otherwise Poisson 
model is thought to take the form of a gamma shape or dis¬ 
tribution, i.e., a(rf>) = 1/k, with k > 0. The NB probability 
distribution function is then given by: 


f(y;k,u) 


r(y + k) (_k_\ h ( 1 _J L _y (7) 

F(k)T(y + l) \p + kj \ p + kj 


The distribution function has two parameters, p and k, al¬ 
lowing more flexible models than the Poisson distribution. 
The symbol T represents the gamma functioif] The mean 
and variance are given by 

2 

Mean = /x; Variance = p + — = p + ap 2 . (8) 

k 

The NB distribution has distributional assumptions similar 
to the Poisson distribution with the exception that it has a 
dispersion parameter a = 1/k to accommodate wider count 
distribution shapes than allowed by the Poisson model. As 
the dispersion parameter, a, approaches 0, lim ap" — 0 or 

OL— >-0 

lim p 2 /k = 0, then the variance equals the mean which 

k—too 

recovers the Poisson distribution. 

It should be noted that if different clusters of counts 
have different gamma shapes, indicating differing degrees of 
correlation within data, and if the NB Pearson \ 2 dispersion 
statistic is greater than one, then the NB model may itself 
be overdispersed; i.e the data may be both Poisson and NB 
overdispersed. Random effects and mixed effects Poisson and 
NB models are then reasonable alternatives (Hilbe|2014|. 

An additional situation should also be mentioned. If the 
Poisson dispersion statistic is less than one, this is evidence 
of Poisson under-dispersed data. The NB model is not ap¬ 
propriate for handling Poisson under-dispersion; however, 
the generalized Poisson model is. We do not discuss under¬ 
dispersed data in this article, but the subject warrants future 
study as to how it applies to astrophysical data. To guide the 
reader, Table [l] displays the main assumptions of the OLS, 
OLS with a log-transformed response variable, Poisson, and 
NB regression models discussed in the previous sections. 


4 DATASET 


As a study case, we use the c atalogue of globular cluste rs 
presented in 


Harris et al. 


120131 (see also 


Harris et al. 


2014 P| 


The data are composed of 422 galaxies with published mea¬ 
surements of their globular cluster populations. There is a 


Table 2. Summary of the parameters used in this work from the 
catalogue of globular clusters compiled by|Harris et al.| 


Parameter 

Definition 

Ngc 

Number of globular clusters 

My 

Absolute visual magnitude 

a 

Bulge velocity dispersion 

M B h 

Central black hole mass 

Mdyn 

Dynamical mass 

€ Nqc 

Uncertainty in Nqq 

e M v 

Uncertainty in My 

€cr 

Uncertainty in a 

€m bh 

Uncertainty in Mbh 


range of galaxy morphologies from which we indexed 247 as 
elliptical (E), 94 as lenticular (SO), 55 as spirals (S) and 26 as 
irregulars (Irr) galaxies for illustrative purposes. Note that 
the original catalogue presents 69 different subcategories of 
morphological classifications which will be discussed in sec¬ 
tion [6] This is a compilation of literature data from a vari¬ 
ety of sources obtained with the Hubble Space Telescope as 
well as a wide range of other ground based facilities. Beyond 
Ngc, we select the following properties for our analysis: cen¬ 
tral black hole mass, dynamical bulge mass, bulge velocity 
dispersion, and absolute visual magnitude as described in 
Tabled 


5 MODELLING THE POPULATION SIZE OF 

GLOBULAR CLUSTERS 

Within this section we demonstrate the application of 
Bayesian GLM regression for modelling the relationship be¬ 
tween Ngc and the following galaxy properties: Mbh, a. My 
and Md yn . Hereafter, unless otherwise stated, the analysis is 
made using a sub-sample of 45 objects from which we have 
observations for all the property predictors. In section [A4| an 
additional analysis uses the entirety of the available data. 

A few common terms in statistical modelling need to 
be reviewed to facilitate our model applications. The analy¬ 
sis focus is the prediction of Ngc as a function of the global 
galaxy properties. Therefore, Ngc represents the count (i.e., 
a non-negative integer) response variable, while My, Mbh 
and Mdyn are interchangeably called covariates, explana¬ 
tory variables or predictors. If included in the model, the 
galaxy morphological type is also considered a nominal cat¬ 
egorical predictor (see section [3]). The whole analysis is per¬ 
formed using JAGS (Just Another Gibbs Sampler)]^] a pro¬ 
gram for analysis of Bayesian hierarchical models using a 
Markov Chain Monte Carlo (MCMC) frameworl0 For each 


3 If n is a positive integer, T(n) = (n — 1)!. 

4 The complete catalogue can be obtained at http://www. 
physics.mcmaster.ca/~harris/GCS_table.txt 


° http://CRAN.R-project.org/package=rjags 
6 Note that count models can be approached by other methods, 
such as a full maximum likelihood algorithm (see Hilbe! 2011, for 
a review) 
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Type +E^S + SO 



Figure 1 . Total number of globular clusters, Nqq, plotted versus 
the central black hole mass, Mbh- The dashed line represents 
the expected value of Nqq for each value of Mbh using Poisson 
GLM regression, while the shaded areas depicts 50%, 95%, and 
99% prediction intervals. Galaxy types are coded by shape and 
colour as follows: Ellipticals (E; blue solid circles), spirals (S; red 
open triangles), and lenticulars (SO; orange asterisks). An ArcSinh 
transformation is applied in the y-axis for better visualization of 
the whole range of Nqq values, including the null ones. 


(“l.®l) ( a 2’) 



NGC;i 


Figure 2. A graphical model of equation © representing the 
hierarchy of dependencies for a data set of galaxies indexed by 
the subscript i. The sinusoidal curves represent stochastic depen¬ 
dencies, while straight arrows a deterministic ones. To save space, 
we replace Mbh by M . in the diagram. 


regression case, we initiate three Markov chains by starting 
the Gibbs samples at different initial values sampled from a 
normal distribution with zero mean and standard deviation 
of 10. The initial adapting and burning phases were set to 
22,000 steps followed subsequently by 50,000 steps, which 
was sufficient to guarantee convergence of each chain for all 
studied cases. 

We now use the relationship between Mbh and Nqc as 
an example to illustrate how the statistical model is built. 
To motivate the use of the more general NB distribution, 
we start the analysis assuming a GLM Poisson regression 
model neglecting the uncertainties in measurements at this 
stage for simplicitjQ This leads to the following model: 


— 

NGC-,i ~ Poisson(/ii); 

Mi = e Vi ; 

r/i — A) + Pi X M BH ;i; 

/3o ~ N(0, 10 6 ); 

Pi~N( 0,10 6 ); 

* = !,■•■ ,N. 


(9) 


This set of equations reads as follows: each galaxy in the 
dataset, composed of N objects, has its globular cluster pop¬ 
ulation sampled from a Poisson distribution whose expected 
value, M) relates to the central black hole mass through a 
linear relation expressed by V- Since we don’t have previ¬ 
ous information about the values of the coefficients /3o and 
/3i, we assigned non-informative Gaussian priors with zero 
mean and standard deviation equal to 10®. We refer the 
reader to appendix [A] for an example of how to implement 
a Poisson GLM in JAGS. The fitted curve for this model is 
displayed in Fig. [l] The grey shaded areas represent 50%, 
95%, and 99% prediction intervals, which are the regions 
where a future observation will fall with these given proba¬ 
bilities Note that the areas in the plot are too narrow to be 
visually discriminated. A visual inspection clearly indicates 
that the Poisson model isn’t adequate to explain the data 
variability since most of the data fall outside the three pre¬ 
diction intervals. Also, the dispersion statistic for this model 
is T) = 1039, which is a strong indication of an inadequate 
model. All other covariates, a, and Mv and Md yn , lead to 
models with similarly high levels of Poisson overdispersion. 
Hence, hereafter we discuss construction of the full model 
based on the NB family to mitigate overdispersion and to 
include the uncertainties in the observational quantities. Un¬ 
like the Poisson model, by employing a NB distribution we 
allow the incidence rate of globular clusters to be itself a 
random variable. 

Continuing with our working example, we keep the dis¬ 
cussion using the relationship between Ngc and Mbh, but 
see appendix [B] for descriptions of the other models. The 
first step is to understand how to include information about 


7 Neglecting the errors at this point does not affect the conclu¬ 
sions regarding the level of Poisson overdispersion. 

8 Not to be confused with the commonly used confidence interval 
in frequentist statistics. A 95% confidence interval will contain 
the sample mean with 95% probability. In other words, a larger 
number of repeated samples from the data would contain the 
sample mean 95% of the time. 
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the uncertainties in the measurements (see e.g., Andreon & 
|Hurn|[2013| for a review of measurement errors in astron¬ 
omy) . Measurement errors in the response count variable are 
the trickiest part to be modelled. The classical model with 
an additive error term y = y* ± e is inappropriate since it 
does not ensure that the observed value y is non-negative. 
The appropriate model is described below and its graphical 
representations are displayed in Fig. [2] 



( 10 ) 


The above is slightly more complex than the model displayed 
in equation (|9| and reads as follows. Each galaxy in the 
dataset with N objects, has its globular cluster population 
sampled from a NB distribution whose expected value, y, 
relates to the central black hole mass through the linear 
predictor y. The additional transformation pi = k/(k + fj,i) 
is required due to how the NB distribution is parametrized 
in JAGS. The uncertainties related to the counts, ejv Gc; i, are 
taken to be associated with the mean, fi, of the NB distribu¬ 
tion and are modelled using a shifted binomial distribution, 
B, with zero mean and taking on integer values in the range 
[—ejv Gc; i,+ejv Gc; i] (see e.g., Chapter 13 from Cameron & 


Trivcdi 2013, from which this approach is loosely based.). 
Uncertainties associated with the observed predictor 
are modelled using a Gaussian distribution with unobserved 
mean given by the “true black hole mass”, Mg H . t , and stan¬ 
dard deviations given by the reported uncertainties in the 
observed black hole mass, &m bh ;i ■ Since Mgg.i is itself an 
unobserved variable, we add a non-informative F prior on 
top of which we added non-informative hyperpriors for the 
shape, ao, and rate, 9o, parameters of the F distribution. 
The choice of a T prior is motivated by the fact that the 
black hole mass is a continuous, but non-negative quantity 
which makes T a more suitable distribution. For the shape 
parameter k, we assigned a non-informative uniform prior, 
W, as suggested in Zuur et al.| ([2013). For the coefficients 
/3o and /?i we assigned non-informative Gaussian priors with 
zero mean and standard deviation equal to 10®. 

Adapting the model above for each combination of Ngc 
and a given galaxy property generates the fitted curves dis¬ 
played in Fig. [3] The grey shaded area represents 50%, 95%, 
and 99% prediction intervals, while the dashed line repre¬ 
sents the expected value of Ngc for each value of the co¬ 
variate. Note the remarkable agreement between the model 


and the observed values with prediction intervals enclosing 
the entirety of the data, including objects that have been 
previously declared outliers and even removed from analy¬ 


sis, such as our own Milky Way (e.g., Burkert & Tremaine 


|2010| |Harris fe Harris|2011| |Harris et al.|2014[). 


5.1 Fit diagnostics 


If the Markov chains are all representative of the posterior 
distribution of the fitted parameters, they should overlap 
each other. Traceplots, Fig. [4] and density plots, Fig. [5] are 
two useful visual diagnostics that are commonly used to test 
for chain convergence. We can see that the chains do mix well 
after the burn-in period, suggesting that the chains are pro¬ 
ducing representative values from the posterior distribution 
for /3o, /3i and k. Additionally, we used a more quantitative 
check, viz., the so-called Gelman-Rubin statistic ( |Gelman] 
& Rubin 19921. The underlying idea is that if the chains 


have reached convergence, the average difference between 
the chains should be similar to the average difference across 
steps within the chains. The statistic equals unity if the 
chains are fully converged. As a rule, values above 1.1 in¬ 
dicate that the chains have failed to properly converge. The 
Gelman-Rubin statistic fell below 1.05 for all estimated pa¬ 
rameters in our analysis. Hence, once we convince ourselves 
that the model is working properly, the next step in the 
analysis is to add interpretations to the fitted coefficients as 
we discuss now. 


5.2 Interpretation of the coefficients 

The exponentiated coefficients e^' of Poisson and NB regres¬ 
sions are also known as rate ratios, or incidence rate ratios, 
which quantify how an increase of unity in the predictor 
variable affects the number of occurrences of the response 
variable. From Table [3] displaying the means and respec¬ 
tive 95% credible intervals of the posterior distribution for 
each parameter, the exponentiated coefficient /3i = 1.59 of 
the Mbh predictor gives a rate ratio of e 1 ' 59 = 4.9. There¬ 
fore, according to the model, a galaxy whose central black 
hole has a mass of, e.g., sa 10 8 Mq has on average approxi¬ 
mately five times more globular clusters than a galaxy whose 
Mbh ~ IO'AFqI^] In other words, one de^°] variation in¬ 
crease in the Mbh leads to an approximately five times in¬ 
crease in the incidence of globular clusters in a given galaxy. 
Likewise, an increase of one dex in Md yn leads to an increase 
of e 219 = 8.9 times in the population size of globular clus¬ 
ters. Another way to state this is, given two galaxies with a 
difference in dynamical mass of one dex, the more massive 
one has a production rate of globular clusters 8.9 times more 
efficient on average. Similar interpretation can be made on 
the other parameters. Another question of interest is how to 
determine the best predictor of Ngc- In the following, we 
discuss how to address this problem from a Bayesian per¬ 
spective. 


9 Note that the analysis was made using log Mbh- 

10 A dex difference of a given quantity x is a change by a factor 
of 10*. 
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Type ~E = S-SO 


Type *E^S^SO 


10 5 --—— 10 5 



Type « E A S^SO Type*-E^S*SO 



Figure 3. Globular cluster population, Nqq plotted against visual absolute magnitude ( My ; top left panel), black hole mass (Mbh 5 
top right panel), dynamical mass (Md yn ; bottom left panel), and bulge velocity dispersion (cr; bottom right panel). In each panel the 
dashed line represents the expected value of Nqq for each value of the covariate using negative binomial GLM regression, while the 
shaded areas depicts 50%, 95% , and 99% prediction intervals. Galaxy types are coded by shape and colour as follows: Ellipticals (E; blue 
solid circles), spirals (S; red open triangles), and lenticulars (SO; orange asterisks). An ArcSinh transformation is applied in the y-axis 
for better visualization of the whole range of Nqc values, including the null ones. 


Table 3. /3i coefficients and scale parameter, k, from Bayesian 
negative binomial regression analysis with Nqc as the response 
variable and Mbh, A7)j yll , cr and My as predictors. The upper 
and lower limits encloses 95% of the credible intervals around the 
posterior means. 


Predictor 

00 

01 

k 


-6.49 ±2.6 

1.59 ±0.30 

1.53 ±0.60 

Mdyn 

-17.72 ±2.75 

2.19 ±0.24 

2.46 ±0.97 

a 

2.99 ±0.78 

0.02 ±0.003 

1.52 ±0.59 

My 

-20.50 ±3.9 

-1.28 ±0.17 

2.23 ± 1.1 


5.3 Model Comparison 

To find the best predictors for the globular cluster popula¬ 
tion, we compare the models using the dispersion statistics 
T> defined i n section|3.2[ and th e deviance information crite¬ 
rion (DIC; Spicgclh alter et al.||2002| ). The latter represents 
a compromise between the goodness of fit and model com¬ 
plexity. It is defined as: 

DIC = Dev + pD , (11) 

where the Dev is the average of the deviance Dev(9) defined 
as Dev(9) = — 21og£(data|6 l ), with L representing the like¬ 
lihood function. The effective number of parameters, Pd, is 
calculated as: 

Pd = Dev — Dev(9), (12) 
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42000 52000 62000 72000 

Iteration 


Figure 4. Illustration of MCMC diagnostics. Three chains were 
generated by starting the Gibbs algorithm at different initial val¬ 
ues sampled from a normal distribution with zero mean and stan¬ 
dard deviation 10. Steps 42,000-72,000 are shown here. The figure 
displays the results for the model .Yqc vs Mbh, with the trace- 
plots for 0o, 01 and k displayed from top to bottom. 


where 9 is the vector of model parameters (/3o,l3i,k for the 
case in study here). The preferred model has the smallest 
value for the DIC statistic. Figs. [6] and [7] depict the results 
for the model comparison using the same dataset. The black 
hole mass displays the lowest values for V and DIC, with dis¬ 
persion statistics as low as T> = 1.05. Although derived from 
an independent analysis, these findings corroborate previ¬ 
ous claims about the tight connection between the central 


black hole mass and globular cluster population (Burkert 


|fe Tre maine 2010|. Nevertheless, it’s worth noting that this 
is not in agreement with a previous analysis performed by 


Harris et al. (2013) using the same catalogue, where they 
found Mdyn as a better predictor for Nqc than MbiRI 


5.4 Further analysis with the entire data set 

Hereafter, we provide a more extensive analysis using the 
entire catalogue of 422 galaxies. The only quantities avail¬ 
able for all objects are the Ngc , galaxy morphological type 
and Mv (H arris et al.| [2013). The advantage of using count 
models for this type of analysis is apparent from the six 
galaxies for which no globular clusters were detected. Such 


ft) 



1.25 1.50 1.75 2.00 

k 



Parameter value 

Figure 5. Overlapped density plots with different colors by chain. 
The plot is a comparison of the target distribution by each chain, 
representing a visual test for convergence. The figure displays the 
results for the model Nqq vs Mbh 5 with the posteriors for /3q, 
and k displayed from top to bottom. 


o 


M v 


Mdyn — 


Mbh 


07 


0.8 



I 


(T9 i tl 1 .2 1 .3 

Dispersion Statistics 


Figure 6. Dispersion statistics, T >, for each model. Values above 
1 represent overdispersion, while values below 1 indicate under¬ 
dispersion. 


11 It is important to note that we are not modelling the same 
relationship as |Harris et al.| who modelled log NgCi the logarithm 
transformation of TVqc , while we model Nqq in the original scale. 
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Figure 7. Deviance information criterion, DIC, for each model. 
Smaller DIC values correspond to preferred models. 

a scenario is naturally accommodated by discrete likelihoods 
while avoiding the failings of logarithmic transformations to 
the response. The statistical model we use is the same as 
that discussed in the beginning of this section and can be 
described as: 



Overall, the model is similar to the one described in equa¬ 
tion (TO). The difference is in the prior for the unobserved 
true absolute visual magnitude, My :i , to which we assigned 
a uniform prior over the range of magnitudes covered by 
the catalogue. The fitted model shows remarkable agree¬ 
ment with the data as displayed in Fig. [8] Very few ob¬ 
jects fall outside the prediction intervals over a wide range of 
galaxy brightnesses. The dispersion statistics for this model 
is T> = 1.15, and the credible intervals for the fitted /3 co¬ 
efficients and scaling parameter, k, are shown in Fig. [9] 
Likewise, as in the previous section, we can interpret the 
f3 coefficient as follows. The mean value of /3i exponenti¬ 
ated is « 0.4, which implies that a galaxy whose absolute 
visual magnitude is one unit greater than another reference 
galaxy has on average 0.4 times less globular clusters; i.e., a 
galaxy brighter by one magnitude over another has on aver- 


Type »E*S S0-»lrr 


10 5 



M v 


Figure 8. Globular cluster population, .Yqc plotted against vi¬ 
sual absolute magnitude My. The dashed line represents the ex¬ 
pected value of A r CC f° r each value of My, while the shaded areas 
depicts 50%, 95%, and 99% prediction intervals. Galaxy types are 
coded by shape and colour as follows: Ellipticals (E; blue solid 
circles), spirals (S; red open triangles), lenticulars (SO; orange 
asterisks), and irregulars (Irr; green open circles). An ArcSinh 
transformation is applied in the y-axis for better visualization of 
the whole range of Nqq values, including the null ones. 


age 2.5 times more globular clusters. Likewise, a galaxy with 
My = —20 has on average 2.5 s « 100 times more globular 
clusters than a galaxy with My = —15, which is consistent 
with a visual inspection of Fig. [8] Another advantage of our 
approach is the possibility to extrapolate the regression so¬ 
lution without making non-physical predictions. The fitted 
model predicts a nearly zero occurrence of globular clusters 
for galaxies with My ^ —11. Considering the total galaxy 
luminosity, L = 10°' 4(Ml O M ' ^Lq, with My Q = 4.83, the 
model suggests that galaxies with 1 ^ 2 x 10 6 Lq are un¬ 
likely to host populations of globular clusters, thus agree¬ 
ing with the literature (e.g., Ha rris et al.| [2013). The use of 
Bayesian prediction intervals allow us to make some inter¬ 
esting predictions: for instance from Fig. [8j we can state that 
galaxies with luminosities L ^ 8.5x10' Lq (or My Q ^ —15) 
should not contain more than 10 globular clusters with 99% 
probability. 

The analysis performed so far did not account for in¬ 
formation regarding different galaxy morphological types. 
Therefore, we are implicitly assuming a pooled estimate 
(e.g., |Gelman fe Hill||2007[ ): all different galaxy types are 
sampled from the same common distribution ignoring any 
possible variation among them. On the other extreme, per¬ 
forming an independent analysis for each class would mean 
making the assumption that each morphological type is sam¬ 
pled from independent distributions and that variations be¬ 
tween them cannot be combined. In the next section we dis¬ 
cuss a more flexible approach together with a brief overview 
of generalized linear mixed models. 
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Parameter value 

Figure 9. Overlapped density plots with different colors by chain. 
The plot is a comparison of the target distribution by each chain, 
representing a visual test for convergence. The figure displays the 
results for the model JVgc vs T/ y , with the posteriors for /3o, /?i 
and k displayed from top to bottom. 


6 GENERALIZED LINEAR MIXED MODELS 

As our final analysis, we introduce one of the most important 
extensions of the GLM methodology known as generalized 
linear mixed models (GLMMs). In particular, we focus on 
one of the simplest GLMM incarnations known as the ran¬ 
dom intercepts model. The random intercepts model, in our 
context, includes an additional term Q to account for a class 
(galaxy type) specific deviation from the common intercept 

Po- 


T/ij — Po + Pi x + Cj, (14) 

where the index j runs from 1 to 69 representing each of 
the different galaxy subtypes reported in |Harris et al.[ A 
standard approach to modelling Q in a standard linear 
mixed regression model is to assume the conditional nor¬ 
mality of the random intercepts with Q ~ A/”(0,1/r), and 
r ~ r(0.01, 0.01). Our intention in incorporating this ex¬ 
tra term into the model is not to simply adjust the data, 
but rather the aim is to identify any particular galaxy sub- 
type which deviates from the overall population mean. For 
this purpose, we employed a popular method for variable 
selection from a Bayesian perspective known as least ab¬ 
solute shrinkage and selection operator (LASSO) which is 
discussed in the following section. 


6.1 Bayesian LASSO 

The original LASSO regression was proposed by |Tibshirani| 
(1996} to automatically select a relevant subset of predic¬ 
tors in a regression problem by shrinking some coefficients 
towards zero (see also |Uemura et ah] [20T5| for a recent ap¬ 
plication of LASSO for modelling Type la supernovae light 
curves). For a typical linear regression problem: 

Vi = Po + Pixi H-h PpX p + e, (15) 


with e denoting Gaussian noise, LASSO estimates linear re¬ 
gression coefficients P = Po + Pixi + ■ ■ • + P p x p by imposing 
a Z/i-norm penalty in the form: 

{ N ( p \ 2 p 1 

V (j/; - V J ( 16 ) 


where re ^ 0 is a tunable constant that controls the level of 
sparseness of the solution. The number of zero coefficients 
thereby increases as re increase. |Tibshiram| also noted that 
the LASSO estimate has a Bayesian counterpart when the p 
coefficients have a double-exponential prior (i.e., a Laplace 
prior) distribution, 

= 27 exp ( - t0 ’ (17) 


where r = 1 /re. The idea was further developed and is known 
as Bayesian LASSO (see e.g., |Park et~a l.] |2008| . Hereafter, 
we use the LASSO formulation for a slightly different pur¬ 
pose, viz., variable selection for random intercept models 
(see e.g.,[Bernardo et al.[20 11 pg. 165). The underlying idea 
is to discriminate between galaxy types that follow the over¬ 
all population mean, i.e. = 0, and galaxies that require an 
additional adjustment in the intercept, i.e. ^ 0. In order 
to include this information, we replace the linear predictor 
ri by equation (141 and add the following equations in the 
model described by equation 


-\ 

C,j ~ Laplace (0, r ); 

r = 1 /re; 

K~r( 0 . 01 , 0 . 01 ); 

3 = !,”• ,69. 


(18) 


The role of the Laplace prior is to assign more weight to 
regions either near to zero or in the distribution tails as 
compared to a normal prior. A visual inspection on Fig. [lO] 
confirms this notion. For the parameter re, we assigned a 
diffuse (non-informative) gamma hyperprior in the form re ~ 
F(0.01, 0.01), which avoids the need of an ad hoc choice of 
re. Note that other possibilities exist such as, e.g., iteratively 
finding re via cross-validation to maximize predictive power. 

Analysis results are displayed on Fig. ED Overall, it sug¬ 
gests that we do not need to add an additional intercept for 
predicting Agc from My. This is consistent with the fact 
that prediction intervals in Fig. [8] enclose ~ 98.8% of the 
data set without any need of a random intercept. Neverthe¬ 
less, the following galaxy types require systematic adjust¬ 
ments: spirals galaxies with moderate size of nuclear bulge 
(Sb), barred lenticulars (SB0), lenticulars (SO) and dwarf el¬ 
liptical galaxies (dEON and dEIN). MG represents one single 
object. Also, UGC 3274 is the brightest galaxy of the galaxy 
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x 

Figure 10. Illustrative comparison between Laplace and Gaus¬ 
sian priors. The Gaussian distribution is represented by dashed 
lines, while the Laplace distribution by solid lines. For all curves 
we assign a zero mean, and the scale (or standard deviation, cr, 
for the Gaussian case) parameters 0.25 (dark blue lines) and 0.5 
(cyan lines). 


cluster ACO 539 ( [Lin &; Moh r |2004| ). Fig. |12| shows that the 
dEON and dEIN objects have a large number of GCs on 
average when compared to other galaxy types with similar 
luminosities, while the lenticulars have systematically fewer 
GCs than expected for the overall galaxy population. This 
can be quantified by looking at the mean value of £ in Fig. 0 
For SO galaxies the mean value of £ is -0.42 indicating that, 
on average, SO galaxies have 34% (1 — e -0 ' 42 ) fewer GCs 
than other galaxy types in the same range of luminosities. 
Generally speaking, galaxy types with 95% credible intervals 
falling on the right side of the dashed grey vertical line in 
ng. mi have more GCs than the overall population mean, 
while galaxy types on the left side have fewer GCs than 
the population mean. While a detailed investigation of the 
causes of this behaviour is beyond the scope of this work, 
it is important to stress the ability of hierarchical Bayesian 
models to explore the multilevel statistical properties of the 
objects under study in an unified way. 


7 CONCLUSIONS 

We employed a Bayesian negative binomial regression model 
to analyse the population size of globular clusters in the 
presence of galactic attributes such as central black hole 
mass, brightness, and morphological type. Hence, demon¬ 
strating how generalized linear models designed to represent 
count data provide reliable outcomes and interpretations. 
The main scientific results and features of our analysis can 
be summarized as follows: 

• The population size of GC is on average 35% lower on 
SO galaxies if compared to other galaxies with similar lumi¬ 
nosities. 

• The relationship between the number of globular clus- 
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Figure 11. Caterpillar plot for the random intercepts G ver¬ 
sus the subcategories of galaxy morphological classifications. The 
thick and thin horizontal lines represent 90% and 95% credible 
intervals respectively. 


ters and other galaxy properties has more variation than 
expected by a Poisson process, but can be well modelled by 
a negative binomial GLM. 

• The Bayesian modelling herein employed naturally ac¬ 
counts for heteroscedasticity, intrinsic scatter, and errors in 
measurements in both axes (either discrete or continuous). 

• Predicted intervals around the trend for expected Ngc 
envelope the data, including the Milky Way, which was pre¬ 
viously considered an outlier. 

• The random intercepts model (with a Bayesian LASSO) 
applied to the correlation between GC population and 
brightness allow us to account for the presence of 69 dif¬ 
ferent galaxy subcategories of morphological classifications, 
and automatically identifies particular types not following 
the overall population mean. Galaxy types dEIN, dEON, 
E/cD, SO, SbO and Sb show significant deviations from the 
general trend. Based on the sample studied here, we advise 
these types to be further scrutinized in order to clarify if 
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Figure 12. Globular cluster population, Nqq plotted against 
visual absolute magnitude My. The dashed line represents the 
expected value of Nqq f° r each value of My, while the shaded 
areas depicts 50%, 95%, and 99% prediction intervals. Galaxy 
types are coded by colours as follows: Ellipticals (E; blue), spirals 
(S; red), lenticulars (SO; orange), and irregulars (Irr; green). As¬ 
terisks represent galaxies belonging to sub-types whose random 
intercept £ is consistent with zero, while circles represent the ones 
with £ 0. An ArcSinh transformation is applied in the y-axis for 

better visualization of the whole range of Nqq values, including 
the null ones. 


there is any physical mechanism behind such deviations or 
merely an observational biaf^| 

• By employing a hierarchical Bayesian model for the ran¬ 
dom intercepts and unobserved covariates (e.g., true black 
hole mass), we allow the model to borrow strength across 
units. This happens via their joint influence on the poste¬ 
rior estimates of the unknown hyper-parameters. 

• If extrapolated, the fitted model predicts a suppression 
in the presence of GCs for galaxies with luminosities L < 
2 x 10 6 Lq. 


• The central black holes mass is in fact a good predictor 
of the number of GCs. One dex increase in Mbh leads to 
an approximate 5 times increase in the incidence of globular 
clusters. The origin of such correlation it is still a matter 
of debate. One possible explanation is that both properties 
are associated with a common event such as major mergers, 
thus galaxies experimenting a recent major merger should 


have a large Mbh mass and GC populations (e.g., Jahnke & 
|Maccio|20 11). The total mass of GCs and the central black 
hole mass can also correlate with the bulge binding energy 
in elliptical galaxies (e.g., Snyder et al.||20TT{ Saxton et al. 


|2014[ ). Rapid growth of the nuclear black hole of a galaxy 
might be fuelled by a massive inflow of cold gas towards the 
centre of the galaxy. The gas inflow would trigger star forma¬ 
tion and the formation of GCs. Hence, leading to an indirect 
correlation between the total number of GCs and the Mbh- 


Scrutinizing which one among these and other possibilities, 
if any, are responsible for this correlation (causal or not) is 
beyond the purposes of this work. However, it does provide 
a clear example on how the adoption of modern statistical 
methods can point to intriguing astrophysical questions. 

A statistical model is based on an appropriate proba¬ 
bility distribution function assumed to generate or describe 
a data set. Hence, the parameter estimating likelihood func¬ 
tion must specify a probability distribution on the appro¬ 
priate scale under study. Discrete data, and count data in 
particular, are not continuous as are data described by the 
Gaussian distribution. The most appropriate way to model 
count data is by using a discrete probability distribution, 
e.g., a Poisson or negative binomial likelihood, otherwise the 
model will likely be biased and misspecified — the price to 
be paid for employing the wrong likelihood estimator for the 
data of interest. 

Generalized linear models are a cornerstone of modern 
statistics, and an invaluable instrument for astronomical in¬ 
vestigations given their potential application to a variety 
of astronomical problems beyond Gaussian assumptions. A 
prompt integration of these methods into astronomical anal¬ 
yses will allow contemporary statistical techniques to be¬ 
come common practice in the research of 21 st century as¬ 
tronomy. 
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APPENDIX A: JAGS MODEL 

Poisson GLM The basic JAGS syntax for a Poisson GLM 
model: 


1 GLM.pois<-model{ 

2 ^Priors for regression coefficients 

3 

4 beta.0"dnorm(0,0.000001) 

5 beta.1~dnorm( 0 ,0.000001) 

6 

7 ^Poisson GLM Likelihood 

8 

9 for (i in 1: N) { 

10 eta [ i ] <-beta.0+beta . 1 *x [ i ] 
n log(mu[i])<-eta[i] 

12 y [ i ] ~ dpois ( mu [ i ]) 


Negative Binomial GLM The basic JAGS syntax for a 
NB GLM model: 


1 GLM.NBC-model{ 

2 ^Priors for regression coefficients 

3 

4 beta.0"dnorm(0 ,0.000001) 

5 beta. 1~ dnorm( 0,0.000001) 

6 k~ dunif( 0.001,10) 

7 

8 GLM Likelihood 

9 

io for (i in 1:N){ 
n eta [ i ] <-beta.0+beta . 1 *x [ i ] 

12 log (mu [ i]) <-eta [ i] 

13 p[i] <-k/ (k+mu[i]) 

14 y[i]~dnegbin(p[i],k) 


Another approach to fit a NB model in JAGS is via a combi¬ 
nation of a Gamma distribution with a Poisson distribution 
in the form (see e.g.,|Marley fe Wand|20lo| |Hilbe|2011[): 


1 GLM.NB<-model{ 

2 ^Priors for regression coefficients 

3 

4 beta.0"dnorm(0,0.000001) 
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5 beta. 1~ dnorm( 0,0.0OOO01) 

6 k~ dunif( 0.001,10) 

7 

s GLM Likelihood 

9 

io for (i in 1: N) { 
n eta[i]<-beta.0+beta.1*x[i] 

12 log (mu [ i ]) <-eta [i] 

13 rateParm [ i ] <-k/mu [ i ] 

14 g [ i ] " dgamma ( k , r at eParm [ i ] ) 

15 y [i] ~ dpois (g[i] ,k) 



APPENDIX B: BAYESIAN MODEL FOR EACH 
COVARIATE 

Dynamical mass versus globular cluster population 

Bayesian NB GLM model for the relationship between Ngc 
and galaxy dynamical mass Mdyn- Since, there is no informa¬ 
tion about the uncertainties in the measurements of Mdyn, 
we neglect this in this particular model. 



Bulge velocity versus globular cluster population 

Bayesian NB GLM model for the relationship between Ngc 
and bulge dispersion velocity a. 
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